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Abstract 

We study metapopulation models for the spread of epidemics in which different sub- 
populations (cities) are connected by fluxes of individuals (travelers). This frame- 
work allows to describe the spread of a disease on a large scale and we focus here 
on the computation of the arrival time of a disease as a function of the properties of 
the seed of the epidemics and of the characteristics of the network connecting the 
various subpopulations. Using analytical and numerical arguments, we introduce 
an easily computable quantity which approximates this average arrival time. We 
show on the example of a disease spread on the world-wide airport network that 
this quantity predicts with a good accuracy the order of arrival of the disease in the 
various subpopulations in each realization of epidemic scenario, and not only for an 
average over realizations. Finally, this quantity might be useful in the identification 
of the dominant paths of the disease spread. 
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1 Introduction 



In our modem world, the existence of various transportation means has strongly 
affected the way in which infectious diseases spread among humans. In fact, it 
has become unavoidable to take into account in the study of the geographical 
spread of epidemics the various long-range heterogeneous connections typical 
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of modern transportation networks, This naturally gives rise to a very com- 
plicat ed evolution of epidemics characterized by heterogeneous o utbreaks pat- 
terns fICohenl . I2OOOI : ICliff fc Haggettl . l2004l : IColizza et aU . hooeM , as recently 



documented at the world-wide level in the SARS case ( ,http:// www.who.int / csr/sars/en ) 



Such complex phenomenon can be tackled at different granularity levels, in- 
cluding very detailed agent-based simulati ons, complicated social and spatial 

structures, complex contact networks, etc. ([Anderson fc Mavl.ll992l:lHethcote &: Yorkd . 



1984j:lKretzschmar fc Morrisl.ll996l:lKeeling 



£l:lMav fc LlovdT I2OO1I : [Ferguson et all 



1999 



Pastor-Satorras &: Vespignani 



2001a 



2OO3I : 



20031 : iMevers et aU . l2005l : IChowell et al 



Eubank et all |2004J ). In particular, a very important class of models in 



modern epidemiology describes the propagation between different intercon- 
nected subpopulations with the help of the so-called metapopulation models. 
In these models, the disease spreads inside each subpopulation (often assumed 
homogeneously mixed) and is transmitted between different subpopulations by 
a "coupling" depending on the model. In the case of human infectious diseases, 
this coupling is caused by the movements of individuals and depends on the 
transportation means relevant at the spatial scale chosen for the description 
of the spread. 



In the context of large scale spread, air-transportation represents a major 
channel of epidemic propagation, as pointed out in the modehng approach to 



global epidemic diffu sion of Rvachev and Longini (IRvachev fc Longinil . Il985 



Baroyan et all 119691 ) and by similar studies on t he behavior of specific out 
breaks such as pandemic influenza, HIV or SARS ( Longini . 19881: Grais al. 



2003; B rownstein et a/.. 12006 



Flahault fc Valleronl . Il992l : iHufnagel et all 12004 



Colizza et a/.l . l2006d Jq. 12007 d ). In this case, the relevant metapopulation model 
considers as subpopulations the inhabitants of the various world cities and the 
air travel of infectious individuals results in the spread of the disease from one 
city to another. 



Taking into account the global airline transportation infrastructure in metapop 
ulation models has recently beco me possible thanks to the availabil i ty an d 



analysis of large scale databases (IBarrat et all 12004 iGuimera et all l2005l ) 



and to the always increasing computer capacities. The word wide airport net- 
work (WAN) is described by a complex weighted graph, in which the air- 
ports are the vertices and the weighted links represent the presence of direct 
flight connections among them (the weights corresponding to the number of 
available seats on each connection). This network, which consists of more 
than 3000 nodes and 18000 links, displays small-world properties as well as 
strongly heterogeneous topology and traffic properti es. For many different dy- 



namical phenomena occurring on complex networks ([Albert Sz Barabasil. 12002 



Dorogovtsev fc Mended . [2003[ : [Pastor-Satorras fc Vespignanj . [2004 : [Boccaletti et al. 
2006[), the presence of such emerging properties has been shown to imply 



the breakdown of the standard results. This is particularly true for epidemi- 



2 



ology (jPastor-Satorras fc Vespignanil . 12001 30; iMav fc Llovdl . l200lh where 
classical results about epidemic thresholds do not apply if the heterogene- 
ity of the network is too large. This therefore calls for a systematic inves- 
tigation of the epidemic spread in the framework of metapopulation models 
defined on complex networks such as the world airport network. Such in- 
vestigations can be carried out in various parallel and complementary direc- 
tions. On the one hand, the availability of large scale datasets has recently 
allowed for the developr nent of tools for intensive computational epidemiology 
(IColizza et l 2006a.fei. l2007Qh. w hich can be used for example for scenario 



evaluations (jColizza et all l2007d ) and risk assessment. On the other hand, 



these computational tools can also be used for the study of fundamental prop- 
erties of these metapopulation models, and to understand the main mecha- 
nisms of the disease spread and the role of the va r ious co mplex properties of 
the transport network. In particular, IColizza et all (12006 J a) have investigated 
how the properties of the global air travel network affect the heterogeneity and 
the predictability of spreading patterns. Very recently moreover, the effect on 
the epidemic threshold of both the network properties and the endogenous 
reaction was studied using arialytical and numerical t ools t ypical of statistical 
physics (IColizza et all l2007a : IColizza &: Vespignanil . 120071 ) . 



In this paper, we tackle the problem of determining the arrival time of an epi- 
demics spread in a metapopulation model consisting in subpopulations (cities) 
coupled by a transportation netw ork. We stress that this study is different from 
the ones of lCrepev efall (l2006h and iBarthelemv et all (l2005f ). in which each 
node of the network is an individual (and not a full subpopulation), and where 
links represent the possibility of the contamination between individuals. The 
main result of the present paper is to exhibit an easily computable quantity 
depending only on the networks characteristics and which approximates the 
average arrival time. This quantity is only an approximation of the exact ar- 
rival time but we show on the example of the world-wide disease spread that 
it allows to rank the cities according to the arrival of the disease. 



The paper is organized as follows. We present in section 2 the Rvachev-Longini 
metapopulation model for the epidemics spreading, and in section 3, we study 
the arrival time of a disease for simple topologies such as a one- dimensional line 
of cities. We then extend these results to the case of complex networks and we 
propose a quantity which is a good approximation to the average arrival time in 
the various nodes or subpopulations connected through an arbitrary network 
of connections. We finally show in section 4, using numerical simulations of 
epidemic spreading on the worldwide airport network, that our approximation 
is not only valid for the average arrival time, but is also able to give the order 
of arrival of the disease in the various cities with a good precision for each 
spreading realiza tion. A brief account of some of the results presented here 



can be found in (iGautreau et all 120071 ). 
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2 The Rvachev-Longini metapopulation model 



The Rvachev-Longini model (IRvachev fc Longinil . 119851 ) was initially intro- 
duced to describe the spread of the 1968-69 Hong-Kong flu at the worldwide 
level. This metapopulation model uses two different levels of description: the 
various subpopulations are given by the inhabitants of the 58 large cities cor- 
responding to the 58 largest airports, and homogeneous mixing is assumed 
at the individual city level. In this model, the propagation of the disease 
from one subpopulation to another is due to in dividuals traveling on the air- 
transportation network between these airports (IRvachev &: Longinil . Il985l ). 



The evolution of the number Jj of infectious individuals in each city i can thus 
be written as the sum of two terms 



dj, = K{{x,}) + n{{i,}) , 



where the flrst term K of the rhs describes the (epidemic) reaction process in- 
side each subpopulation (city), due to the interaction of individuals in the 
various possible states X {X = S,L,I,R... depending on the population 
compartimentalization into Susceptible, Latent, Infected, Recovered... indi- 
viduals). The second term of the rhs represents the incoming and outgoing 
fluxes of infectious individuals to and from other cities j. This model there- 
fore considers a simplifled mechanistic approach with a markovian assumption 
in which individuals are not labeled according to their original subpopulation, 
and where at each time step the same traveling probability applies to all 
indivi duals in the subpopulation, without any memory of their previous lo ca- 
tions ( Rvachev &: Longini . 1985; Hufnagel et al . 2004 : Colizza et al . 2006 d lil). 
The travel term Q depends therefore on the air-transportation network: if the 
weight Wij represents the number of passengers traveling from i to j per unit 
of time {wij = if there is no direct connection between i and j), and Ni the 
population of the city i, it is reasonable to assume that the probability per 
unit time that an individual in city i travels to city j is given by Wij/Ni. In 
the case of a simple SI model, the time evolution of the number Jj of infectious 
inflviduals in city i is therefore given by: 



dJi = Mi (t) 



- I, (t) 



-E 



(2) 



where A is the transmission rate. Similar equations can be written for the other 
compartments in the various possible models (SIS, SIR...), by modifying the 
flrst term K of the rhs in Eq. ([2]) accordingly. 



This original formulation has an important drawback. Indeed, it is determin- 
istic, since it considers only expectation values, while both epidemic spreading 
and travel of individuals are inherently stochastic processes. The number of 
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infectious individuals is then treated as a continuous variable and, even if the 
initial condition (at t = 0) of the spreading consists in one single infectious 
individual in a given city io, all cities have a non-zero density of infectious at 
any time t > 0. In order to tackle the problem of the arrival time of the spread 
in the various subpopulations, t his pathology has to be ac counted for. A first 
possibility, already presented by lRvachev fc Longinil (Il985l ). to avoid this unre- 
alistic situation consists in considering that each travel term Wijli/Ni is present 
only when WijIi/Ni > 1, i.e. in putting a threshold term 6{wijIi/Ni — 1), cor- 
responding to the fact that the expected number of travelers has to be larger 
than unity. In this approach however the "quantity" of individuals traveling is 
still a continuous variable. Another possibility, which consists in treating all 
quantities as integers and going back to the micr oscopic stochastic pro cesses, 
can be necessary for detailed modeling purposes (jColizza et all l2007ol ) . This 
solution is however computationally quite demanding, so that we will consider 
in all numerical simulations performed in this paper an intermediate frame- 
wo rk: the travel t erm is treated as i n the stochastic gener alization described 
by Colizza et al. (200633) (see also Hufnagel et al. ( 20041 )). where the num- 
ber of individuals traveling on each connection is an integer variable randomly 
extracted at each time step. More precisely, each of the [/j] (the integer part 
of Jj) infectious individuals has a probability = WijAt/Ni to go from i 
to j in the time interval At. The numbers of infectious individuals trav- 
eling on the various connections existing from i to other cities j are given 
by a set of stochastic variables which follow then a multinomial distribution 
flColizza et di 12006 d ) 



P ({e.,} I {/.}) 



X T[p 



(3) 



where the mean and variance of the stochastic variables are ([/«])) = 
WijAt[h]/Ni and Var {^ij) = ^At (l - ^At) [/,]. On the other hand, the 
endogenous growth (inside each city) will be treated for simplicity as a de- 
terministic evolution of densities of infectious and susceptible individuals. We 
have moreover checked numerically that the inclus ion of stochast i c effec ts de- 
scribed by noise terms in the evolution equations ( Colizza et al . 2006 d J^ do 
not change our results. Finally, we will present the results essentially for the 
simple SI model, since arrival times will basically depend on the first stage of 
the disease development. We have however carried out numerical simulations 
as well in the SIS and SIR case, with consistent results (see section 4). 



It is interesting to note that in realistic cases such as a disease propagatio n on 
the WAN, most weights are symmetric {wij = w 



. as a disease propagatio n on 
'ji) fiBarrat et ali booi ] but 



the probabilities of travel from one city to another are not {pij ^ pji). The 
travel therefore effectively occurs as a random diffusion on a directed weighted 
complex network. Moreover, an important difference with the case of random 
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walks on complex networks (INoh fc Riegeii . |2003| ) comes from the endogenous 
evolution inside each subpopulation. Indeed, as soon as an infectious individ- 
ual reaches a city, a contamination process starts and the number of random 
infectious walkers is not constant. Such an intricate behavior has also impor- 
ta nt consequences on th e existence of epidemic thresholds, as studied recently 
bv IColizza eToZI t007Jh . 



As a last general remark, we note that the simple topological distance between 
nodes should naturally play a role in the context of arrival times, but does not 
contain all the information needed to characterize such a process, nor does 
a priori the opti mal weighted distance, which takes into account the weights 
( IWu et all I2OO6I ) but not the populations nor the endogenous epidemic evolu- 
tion. Moreover, since most transportation networks are small-world networks 
(even when taking into account weights), many cities lie at the same topolog- 
ical distance from a given seed, but will be reached at very different times. 



3 Arrival times statistics 



In this section, we study the arrival time of the first infectious individual in 
the various cities, starting with very simple topologies for the transportation 
network. 



3.1 Two cities : arrival time distribution 



We start with the case of two connected cities (0 and 1). We denote the link 
weight by wqi = w, the population of the first city by Nq = N, the number 
of infectious individuals in city by Jq = /. The probability to jump from 
to 1 during the time interval At is p = ^At. The initial condition is given 
by /o(t = 0) = /°, Ii{t = 0) = 0, i.e., infectious individuals are present only 
in city 0. We consider that the travel events occur as instantaneous jumps (of 
probability p for each individual) at discretized times. The probability that 
the first infectious individual arrives at time ti = t = nAt in city 1 is then 
given by 

n-l 

P (t, = nAt) = (1 - (1 - p)[^o(«^*)]) X n (1 - p)t'°^^^*^^ , (4) 

i=l 

which expresses the fact that at least one successful "jump" from to 1 of 
an infectious individual occurs at time nAt, and none at previous times. In 
Eq. (jl]), the quantity [Jq] denotes the integer part of Jq. In order to obtain 
the probability density of the arrival time in city 1, we first note that in real- 
world systems the number of travelers is usually small with respect to the 
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total population of a city: p = wAt/N ^ 1 (with any reasonable At). In this 
limit, Eq. (jlj) becomes 

P (ti = t)= p[Io{t)]e~P^o<.<n ^'oi^m _ (5) 

We will also assume that the travel probability is sufficiently small so that the 
number of infectious individuals Io{t) in the city grows substantially before 
the city 1 is contaminated. Writing 1 -C Jo(t) -C (the second inequality is 
due to the fact that in realistic cases only small fractions of the population 
are contaminated), we thus consider the following approximation: [I]{t) ~ 
/o(t) ~ /°e'^* (as long as t < ti). Using as well the standard approximation 
AtY.o<i<n[Io (iAt)] = /o Jo (r) dr, we obtain for J° = 1 

P{t)dt = ^e^'-wx-''e(t)dt , (6) 

where 6(t) is the Heavyside function which ensures the positivity of the ar- 
rival time. The distribution ([6]) is a Gumbel distribution with average (ti) = 
J (\n (^^^ — 7), where 7 is the Euler constant. The contribution of the nega- 
tive t in the distribution has to be negligible before 1, which reads P (t) dt = 
<^ 1. This range of validity corresponds also to the hypothesis of A(ti) ^ 1 
(the time to increase /q by one individual is small with respect to the time 
scale of the epidemic arrival in city 1). 

Figure [T] displays the results of numerical simulations of the propagation be- 
tween two cities, using discrete stochastic travel events as described in section 
2. The whole distribution of arrival times in the second city shows a very good 
agreement with the form ([6]), obtained through the continuous approximations 
detailed above. When w/{NX) is not small enough, stronger deviations are ob- 
tained; this is expected since the hypothesis 1 <^ /o(t) for t < ti is less valid. 
Nonetheless, the overall shape of the distribution is still in good agreement 
(Inset of Fig. p. 



3.2 One- dimensional line 

We now consider a one- dimensional line of subpopulations ^ > 0, with an 
initial condition given by one infectious individual in city 0. We denote by Wi 
the travel flux between city i and city i + 1, tj the arrival time of the ffist 
infectious in city i, and Aj = tj — tj_i (tg = 0). The quantity t„ is thus the 
sum of the random variables Aj for i = 1, ...,n, and the average is given by 
(tn) = J2i=i....,n{^i)- These random variables Aj are a priori correlated and 
not identically distributed which implies that the central limit theorem can 
not be used in general. 
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Homogeneous line 



In the case of a homogeneous hne with uniform populations and weights {Ni = 
N, Wi = w), one expects that the distributions of A„ becomes independent of n 
at large n, with a well-defined average lim„_^oo(An) = {^), so that t„ ~ n{A). 
The issue then is the computation of (A). Two time scales only are present: 
N/w which represents the average time for an individual to travel from one 
city to the next, and 1/A which represents the typical transmission time of 
the disease. Dimensional analysis thus implies that the adimensional quantity 
A (A) has the form 

In order to estimate the unknown function F, a first possibility is to consider 
that at large n the spread consists in an evolving epidemic front obeying the 
continuous equation dtl{x, t) = \SI + wdll{x, t). Looking for a traveling wave 



solution (IMurrayl . 120041 ) leads to a front speed v = 2w and thus to 



I NX 
(AW--W— . 



Another possibility consists in using the results of the previous subsection, 
and assume that the average of A„ remains close to the arrival time in the 
first city (A„) ^ (Ai) = (ti) which yields 



/A\ _ 1 

{A/Gumbel — 




-7 



(9) 



This approximation neglects the fact that increases between and 

t„ due both to the endogeneous growth in city n — 1 and to the arrival of 
infectious individuals from city n — 2, while for the computation of (ti), only 
the endogeneous growth of Iq has to be considered. It can therefore be expected 
that (A)^^^^; will overestimate the real (A). 

Finally, we can also consider a deterministic formulation, in which travel be- 
tween i and i + 1 occurs only if there are enough infectious individuals in i, i.e. 
Ii{t) has reached a certain threshold 6, and is then treated as continuous. In 
this framework, each city is considered infected only if the number of infectious 
individuals is above 6, and tn is defined as the first time that /„ reaches the 
threshold: In{tn) = Between tn-i and t„, no travel can therefore occur out 
of n, and we can write dtin = XIn + ^In-i- During [t„_i,tn], under the real- 
istic hypothesis that /„ <C Nn, and that is in a first approximation given 
by In-i = e^*^*"*"-^'', we obtain the equation 1 = ^A^e'*'^". The difference 
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between the arrival times in two successive cities is therefore given by 




) 



(10) 



where W is known as the Lambert W function. 

In order to test these various analytical approaches, we consider numerical 
simulations of the stochastic model described in section 2, for identical cities 
of population = 10^ located on a one- dimensional line, with uniform travel 
fluxes w between successive cities. The measured average arrival times are, as 
expected, proportional to n at large n (not shown). Figure [2] displays the cor- 
responding slope as A(A) versus ln{XN/w). Various values of A^, A and w with 
the same ratio XN/w yield the same A (A), as predicted from the theoretical 
dimensional analysis ([7]). As shown in Fig. [21 (A) front is in agreement with 
the simulations only at small XN/w: the slowest the travel, the less a spa- 
tial continuous approximation is valid. On the other hand, both {A)det and 
{A)Gumbei display reasonable agreement, and in particular correctly capture 
the increase with ln(AA^/tf) for large XN/w. As expected, {A)Gumbei slightly 
overestimates (A). We also show in the inset of Fig.[2]the correlations between 
the A„'s. These correlations vanish rapidly with the distance along the line. 
Negative correlations can however be observed between A„_i and A„. Such 
phenomenon can be understood as follows. Assume that, in a given spreading 
realization, A„_i is small. In this case, In-2 {tn-i) will be unusually small too. 
The "reservoir" of infectiousness defined by city n — 2 will thus transmit less 
infectious individuals to n — 1 at times t > t„_i, and therefore the subsequent 
contamination of city n will be slower, leading to a larger A„. 

The numerical simulations allow also to measure the whole distributions of ar- 
rival times and of their intervals. In Fig. [31 we show P( A„) for different values 
of n. The Gumbel shape is only valid for n = 1, and for large n more sym- 
metric (Gaussian-like) distributions are obtained. On the other hand. Fig. [H 
shows that the distribution of the arrival times themselves display asymmetric 
Gumbel-like shapes even at large n. 



General populations and weights 

We now consider the general case of heterogeneous populations Ni and fluxes 
Wi- The application of the same considerations as in the case of a homogeneous 
hne leads to (t„) = J2'i=i{'^i)- Since our main goal is to understand the case 
of complex networks which usually have the small-world property (i.e. a small 
diameter varying typically as log A^), we will now focus on the properties of t„ 
at small n. Assuming that the average of Aj remains close to the formula for 
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the arrival time in the first city ((Aj) ^ {ti{X,Wi, Ni))), we obtain 



A(t„) ^ Xn = In 



NiXe-^ 



n-l 

n 

.i=0 



This expression hnks the expected arrival time in city n of a stochastic spread- 
ing event to the quantity Xn which depends on the properties of the line on 
which the spreading propagates. It is symmetric by permutations of weights 
and populations, and Figure [5] shows in fact that the whole distribution of ar- 
rival times, obtained through numerical simulations of the discrete stochastic 
dynamics, and not only its average, respects this symmetry. In the numer- 
ical simulations, heterogeneous populations and travel fluxes are uniformly 
distributed {wi G [10,2000] and Nt e [10^2.10^]). Figure [5] shows that the 
distribution of arrival times is invariant when one replaces (i) all the random 
weights by their geometrical mean w = (Jli=o...nU!iY^'^; (ii) all the random 
populations by their geometrical mean N = (nj=o...n ^i)^''"; (iii) all weights 
by w and all populations by N. The ratios of the average arrival times for 
these different sets {iVj))/(t(w, {iVi})), {t{{wi}, {Ni})) / {t{{wi},N)) 

and {t{{wi}, {Ni})) / {t{w, N)) stay very close to 1, with deviations at most of 
the order of 5%. 

We test in Fig. [6] the relation (fTTl) which predicts the average arrival time 
measured in numerical simulations of the stochastic model. We use for sim- 
plicity uniform populations N = 10®, and we measure, for each fixed set of 
weights, the arrival time in each city, averaged over 10^ stochastic realizations 
of the spread. Figure [6] displays X{ti) versus Xi ^^i the various cities of a 5 
cities line, for 100 different random sets of weights (each point corresponds 
to one given set of weights). We first note that X{ti) is a monotonous func- 
tion of Xi- The figure also clearly shows that x systematically overestimates 
the average arrival time. This is due, as already explained in the discussion 
of Eq. to the fact that the use of {A)Gumbei neglects correlations, and 
in particular the travel between cities n — 2 and n — 1 for the estimate of 
the arrival time in n. Since this travel term increases it increases the 

probability that an infectious travels from n — 1 to n, and therefore reduces 
tn- This overestimation effect is not present for the first group of points that 
lie on the diagonal of the figure and correspond to cities directly connected 
to the seed. For this first group of cities, the average arrival time is indeed 
correctly given by x- On the other hand, the small overestimation can already 
be observed for cities at distance 2 of the seed. For example, points A and 
B highlighted in Fig. [6] have the same Xb = Xa = X- Point A corresponds 
to a city directly connected to the seed S, with a flux of wa passengers per 
unit time, such that xa = Iti{NX/wa) — 7, while point B is obtained for a 
city which lies at topological distance two from the seed, with a city B' in 
between and such that xb = In (^^^77) +ln (^ J^^e-y ^ '^^^ cases correspond 
to the same value of but to different average arrival times: the arrival time 
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in B is smaller than in A, even if xa = Xb, because relation (fTT]) neglects the 
travel of additional infectious individuals between S and B' to compute (ts). 
Despite this systematic effect, Fig. shows that relation (ITTi) allows to define 
a quantity which depends only of the populations and passenger fluxes, and 
such that the average arrival time of the spreading is a monotonous function 
of this quantity. In the following, we will investigate how to generalize this 
quantity to more complex topologies. 



3.3 From the one- dimensional line to complex networks 



Equation flTT]) gives an estimate of the average arrival time in a city n con- 
nected to the seed of the spreading through links with certain fluxes, and 
corresponds to a sum of the quantity ln{NX/w) — 7 along the links followed 
by the disease from the seed to city n. In a network, the sum can be com- 
puted along any of the paths linking the seed s to any other city j. Since these 
different paths correspond a priori to different fluxes and to different arrival 
times, a natural approach consists in approximating the average arrival time 
in j, starting from a given seed s, by the minimum of ffTTl) over all possible 
paths: 

~ = min J2 ' ^''^ 



Inp^ -7 



'12 



where s is the seed, {Psj} is the set of all possible paths connecting s to j, 
and the sum is performed over the links {k, I) along each path. In other terms, 
we have introduced a new weight on each oriented link {k, I) of the network: 
— 7. The quantity x (j k) is then the weighted distance between the 



ln( 

seed s and the node j on this non-symmetric network (such a qua ntity can 



easily be computed for any topology using the Dijkstra algorithm (jPijkstra 



1959h l Note that since the weights are real-valued, it is highly improbable 



that two different paths with the same sum of weights exist, so that equation 
(1121) selects a unique path between s and j. 

Differences between propagation on a one-dimensional line and on a network 
can be expected from the two following important topological differences be- 
tween these structures: (i) paths go through nodes with degree larger than 2 
and (ii) there are usually more than one path between two points. We first fo- 
cus on the effect of the intermediate nodes with possibly large degrees through 
which the disease propagates. To this aim, we consider the simple topology 
consisting of a node connected to k neighbors (as shown in the inset of Fig. [7]). 
The infectious individuals, when they arrive at the central node, have multiple 
possible travel destinations. The probability for each destination to become 
infected is thus decreased; if the seed is in for instance, the initial infectious 
individual may even travel to a different peripherical node before creating 
an endogenous epidemic growth in the hub. This effect can be quantified by 
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measuring the average arrival time of the disease in a given peripheral city i 
in various situations. We first consider the case of a spreading phenomenon 
starting in a randomly chosen peripheral node, 7^ h ci-nd compare the av- 
erage arrival time in i with the average arrival time in the one-dimensional 
situation in which only cities zq, and i are present. The ratio of these two 
times is displayed in Fig. [7] ('seed: not hub') as a function of the degree k of 
the central node. We also show in the same figure the ratio of arrival times 
when the seed is 0, i.e. the hub itself ("seed: hub"), and the topology is either 
the one of the inset, or the simple two-cities configuration in which only and 
i are present and connected. In all cases, the arrival time in i is increased by 
the presence of other possible connections from the hub. The effect is stronger 
for larger degree when the seed is the hub itself. The multiplicity of possi- 
ble destinations there fore not only decreases the predictability of the spread 
(IColizza et all uOOQa ) , but also increases the average disease arrival time in a 
given city. 



Another important point distinguishes a propagation in a network from a one 
dimensional line: namely, that multiple paths can ofte n be found be t ween two 



given nodes. In fact, it has already been shown by ICrepey et al\ (120061 ) in 



the case in which nodes represent individuals and not subpopulations, that 
the average arrival time is decreased by the presence of multiple paths, and 
that not only the shortest paths contribute to the spread between the seed and 
other nodes. To quantify this effect in metapopulation models, we consider the 
simple topology depicted in Fig. [8j a node A (population A^^) is connected 
to a node B by two different paths of length two. The first path connects A 
with weight Wi to an intermediate city C of population A^, which is in its 
turn connected to B with weight W2- The second path has weights w[ and w'2, 
and the intermediate city C has population A^'. We denote the probability 
that the disease spread reaches the city B at time t by P (t) if the two paths 
are present, by Pacb (t) if only the first path through C (weights wi, W2) is 
present, and by Pacb if) if only the second path (through C", weights w'l, 
W2) is present. Since the epidemics reaches B from A through one path or the 
other. Pit) can be expressed as 



P (t) = (1 - ^* Pace (r) d?j Pace (t) + - Pace (r) d?j Pace (t) 

(13) 

The data shown in Fig. H] suggest that we can assume that Pace and Pace 
are Gumbel distributions, of averages Xace = In ('^^) + 1^ ('^) ~ ^'^'^ 
Xace = In ('^^^) {^ut) ~ respectively. We introduce the quantities 

Weq = ^^^l^'' and w'^g = —^^r^- The quantity Weq can be seen as the travel 
flux which, if A and B were directly connected, would yield the same arrival 
time distribution than the real two links AC and CB with weights wi and W2'- 
Xace = In (^) - 7- We then have 
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Weq Ai-^e^' 



Pace (t) = ^6^'-—^'' (14) 



Pac'b (t) = ^e'*-^^ (15) 

and finally from Eq. (fT3l) (using the fact that the weights are small with respect 
to the quantities A^A, A^'A, A^aA) 

P(t) = ^^exp(^At-^^e j . (16) 

We note that P{t) is also a Gumbel distribution, and the average of the arrival 
time in B, (tmp) (mp=multiple paths), can be easily calculated: 

A(tmp) = Xrnp = In ( , ] " 7 , (17) 

SO that the existence of two paths results in the law 

g-Xmp _ g-XACS _|_ g-XAC'S . n 8] 



We have checked this prediction by numerical simulations of stochastic spread- 
ing in the small network formed by nodes A, B, C and C, in the simple case of 
Na = N = N' and wi = W2 = w, w[ =102= w'. Figure [S] displays the ratio of 
the average arrival time in 5, for a spread seeded in A, to the average arrival 
time in B when only one path (the one with weights w) is present. We consider 
only w' < w since the situation with w' > w is obtained by exchanging the 
two paths. As w' —>■ 0, this ratio goes to 1, while, in the other extreme case 
w' —>■ w, a decrease of the average arrival time close to 10% is obtained. More- 
over, the prediction f|T8|) . shown as continuous line, is in excellent agreement 
with the numerical data. 

Interestingly, when more than two paths are present, Eq. f|T8|) can easily be 
generalized to a sum over all possible paths (not necessarily of length 2): 

Mtmp) = Xmp with 

g-Xmp ^ g-Xpath ^ (]^9) 

paths 

where Xpath is computed on each path by Eq. ffTTj) . Figure [9] presents a com- 
parison of numerical results with this analytical prediction, for various path 
multiplicities, showing a very good agreement. Similar results are obtained for 
paths of larger lengths. 

We have also investigated numerically, on simple topologies but with random 
weights, how well the average arrival time in a city is correlated with Xmp 
computed as in (fT9l) . In particular. Fig. [10] displays the average arrival time 
in city B, connected to the seed A of the spreading by k paths of length / and 
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random weights, as a function of both xb given by Eq. (fT2|) (i.e. computed 
along the path that minimizes x) ci^nd Xmp which takes into account all paths. 
While Xmp is slightly more strongly correlated with (t) than xb, and closer 
to its actual value, the improvement is not striking. We have also developed 
an algorithm to compute Xmp on a complex network, by taking into account 
the contributions of all paths of topological lengths d and d + 1 (where d is 
the length of the shortest path, in terms of topology, between the seed and 
the considered city). While the algorithmic complexity is in creased (we used 
techniques based on the Brandes algorithm (Brandes, 200 ll )) with respect to 
the computation of ffT^ . the obtained improvement was again significative but 
not striking. 



4 The world-wide airport network 



We now proceed to numerical simulations of the metapopulation model on 
the Worldwide airp ort networhE] , which displays variou s level s of complexity 



and h eterogeneity (IBarrat et all . l2004l : iGuimera et all l2005l : IColizza et al. 



2006 d l^. At the topological level, the degree distribution is broad and can 



be approximated by a power-law; the links' weights (fluxes) are also broadly 
distributed and span several orders of magnitude. Finally, the world city pop- 
ulations are also broadly distributed according to Zipf's law (jZipj . 119491 ). All 
these heterogeneity levels have been shown to play relevan t roles in the spread 
of epidemics at the worldwide level ( Colizza et al. . 2006 (j lil). We perform the 
numerical simulations of the stochastic spreading on the network containing 
the 2,400 largest airports, which takes into account 98% of the total traffic 
(the WAN comprises 3, 100 airports, but we consider here only the links veri- 
fying In {j^^ > 0, which amounts to the removal of the smallest airports; in 
particular, a certain number of nodes have a larger traffic than inhabitants, 
which leads to a local breakdown of the relation between pij and the local 
population and travel flow). 



4-.1 Average arrival time 



Figure fTTl displays the average arrival time in the various cities, A(tj), as a 
function of Xi as defined by (fT2l) . For each given seed, averages are done over 
1, 000 stochastic realizations of the disease spread. The figure clearly shows 



* We have as well considered artificial networks, with similar numerical results. 
Some arbitrariness in the distribution of weights and city populations is necessary 
for articial networks so that we prefer to present data obtained with a real-world 
network, in which all quantities stem from real data. 
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that the value of x determines the average arrival time (the two quantities are 
very strongly correlated) and various cities with the same x ^-re reached at the 
same time by the propagation. As in the case of the simple topologies studied 
above, x f^ict systematically overestimates the correct average arrival time, 
but can still be considered as a very good approximation. 

Figure [TT] also highlights the effect of hubs, analyzed in the previous section. 
Frankfurt is indeed the node with the largest degree in the network, and the 
arrival times at a given value of x is larger than for other seeds with a smaller 
degree. This effect remains however quite small and is essentially limited to 
the first reached cities. The arrival times seem to be closer to our Ansatz. 
This stems from the delay introduced by the large degree of the initial seed, 
as discussed in the previous section, which increases all the arrival times, 
therefore slightly shifting the data upwards and partially compensating for 
the fact that x overestimates these times. We finally note how very similar 
results are obtained for the SIR model. As expected, the important time scale 
in the SIR (as well as in the SIS) is A— /x: this quantity controls the endogenous 
growth of the epidemic at small times. We have investigated a wide range of 
values for A — /i, as well as for A in the simple SI case. The arrival time is well 
determined by x long as the condition w/NX << 1 given in section 3 is 
satisfied. 



Figure [TT] conveys the result that the average arrival times are not exactly 
given by but are at least determined to a large extent, given a seed, by 
this quantity. An immediate application of this result is given by a possible 
prediction of the order of arrival of the disease in different cities. In order to 
compare the list of cities ranked by the average arrival time and the same 
list ordered according to we compute the indicator known as Kendall's r. 
This indicator allows a quantitative analysis of the correlations between two 



rankings of n objects (INumerical Recipes! . l2004j ) and is given by 



T = ^ 20 

n{n -l)/2 ^ ' 

where Uc is the number of pairs whose order does not change in the two 
different lists and is the number of pairs whose order is inverted. This 
quantity is normalized between —1 and 1: r = 1 corresponds to identical 
rankings while r = is the average for two uncorrelated rankings, and r = — 1 
is a perfect anticorrelation. 

We present in the first column of Table [T] the values of r for the lists of cities 
ordered according respectively to the average arrival time of the spread and 
to the values of x, for various seeds. The second column of the table gives 
the values of r when the lists are ordered by average arrival time and by 
Xmp- In both cases the values are extremely high (two random permutations 
of a list of size N = 2,400 would give an indicator normally distributed with 
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mean and variance g'^i^^^i^ = 2.10"^ (jNumerical Recipea . l2004j )). In order 
to better emphasize this strong correlation between the ordered hst, we also 
display in Table [1] the r of Kendall between the city list ordered by average 
arrival time and by different distances from the propagation seed. As noted in 
the introduction, the topological distance d from the seed is not completely 
irrelevant, but other weighted distances which take into account the diversity 
of fluxes and populations could be expected to be more strongly correlated 
with the arrival time. In particular, a first possibility consists in defining the 
effective "length" of each edge as the inverse of the weight, iij = l/wif the 
disease will spread more easily on an edge if many passengers travel across 
it. The corresponding weighted distance between nodes on the network is 
noted di/w Since the ratios N/w moreover appear naturally in the travel 
probabilities, we also consider that a directed length Ni/wij can be defined 
on each edge, and the distance between each node and the seed can as well 
be computed (dN/w)- The results indicate significant correlations between the 
average arrival times and these various distances. The correlations are however 
much stronger with x and Xmp, showing that these quantities can be used with 
a good confidence as an estimate of the average arrival time order. Very similar 
results are obtained for SIS or SIR models, with respectively Tsis{{t),x} ~ 
0.875 and TsiR{{t), x} ~ 0.866, for a spreading process seeded in Hong-Kong. 



4-2 Arrival order for a given realization 



The previous results are valid for the average arrival times, and it is legitimate 
to ask about their relevance to the case of a single spreading event. Indeed, 
in the real-world there is not such a thing as averages over different realiza- 
tions. The natural extension of our results therefore consists in checking if the 
quantity x can predict the order in which the disease will reach the various 
nodes (cities) of the network in each stochastic realization of the spread. More 
precisely we can compute, for each pair of nodes (i, j) in the network, and 
for a given infection seed, the probability that the arrival times in i and j for 
the same realization of the spread, ti and tj, are correctly ordered by their 
values of i-e. that they verify {ti — tj){xi ~ Xj) > 0. We use the notation 
{h j) =\ Xj ~ Xi \ aiid Fig. [T2] displays the probability fc that a couple of 
nodes with a x-difference Ax{i,j) = Ax are reached by the disease in the 
order predicted by their values of x- If Ax(«,j) = 0, no prediction can be 
made and we obtain indeed fc = 0.5. On the other hand, if Ax{i,j) is large 
(> 10), the two nodes are reached in the correctly predicted order in every 
realization of the spread. 

Since not all node pairs have very different values of x (and thus a large value 
of Ax), we also show as in Fig. [12] the cumulative distribution p>(Ax) of 
the number of couples of nodes with a given value of Ax. For instance, for a 
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spreading process seeded in Hong-Kong (degree 109), 75% of the couples have 
Ax > 2, and these pairs are correctly sorted with a probability larger than 
80%, instead of only 50% on average if no information is available. Figure [T^ 
moreover shows that the precision is higher when the seed has a small degree. 
On the other hand, if the seed is Frankfurt (of degree 240, data not shown), 
the situation is slightly worsened, as can be exp ected from the smalle r global 
predictability of the spread starting from a hub fICohzza et all . l2006(J ) , due to 



the many possible travel destinations available for the infectious individuals. 
In this case, 80% of the couples have Ax > 2, and these pairs are correctly 
sorted with a probability larger than 70%; the probability fc is larger than 
90% for 47% of the pairs. Similar results are obtained when the endogenous 
growth of the disease is described by SIS or SIR models. 



We also present in Fig. [T3] the same quantities fc and p> (Ax) for couples of 
nodes having the same topological distance from the seed. While this topologi- 
cal distance is indeed too simple a quantity to be really useful in the prediction 
of arrival times, it can seem reasonable that a node at distance 1 from the seed 
will anyway be reached before a node at distance 5. If large values of Ax were 
obtained only for nodes at very different topological distances from the seed, 
the prediction shown in Fig. [T2l would not be particularly relevant. We see in 
Fig. [13] that nodes at the same topological distance from the seed can have 
very different values of x and be therefore reached by the disease in an order 
which can be almost certainly predicted by measuring their values of x- For 
example. Fig. [T3] shows that more than 70% of the couples are well ranked 
with a probability larger than 70%, even if the two nodes are at the same 
topological distance {d = 2, 3, 4 or 5) from the seed. Results are presented for 
a seed with large degree, and shghtly better results are obtained for smaller 
degree seeds. 



Finally, another quantitative indication of the relevance of the quantity x is 
obtained, similarly to the previous subsection, by comparing the list of cities 
ordered either by x or by the arrival time of the spread in a given realization 
(and not the average arrival time as in the previous section). Kendall's r is 
therefore now a quantity which fluctuates from one realization to the other (for 
a given seed), and the corresponding histograms are shown, for various seeds, 
in Fig. [TH The values obtained are systematically larger than 0.5, denoting 
a strong correlation between the lists ordered according to the arrival time 
and X- Moreover, the figure highlights how seeds with larger degrees lead to 
smaller values of r. This is in agreement with the slightly better prediction ca- 
pacities displayed in Fig. [T^ when the seed has small degree. Inte restingly, this 



result c an also be related to the issue of predictability studied by lColizza et al. 



(12006 al Jih: the predictability of a spreading process, as measured by the simi- 
larity between two stochastic realizations with the same initial conditions, has 
indeed been shown to be larger when the seed has a small degree. 
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5 Conclusion 



In this paper, we have studied a metapopulation model for the spread of epi- 
demics on a large scale in which subpopulations (cities) are linked by fluxes of 
passengers. Such models are particularly useful for the analysis of epidemics 
which propagate worldwide along the airline conne ct ions (jRvachev fc Longini , 



19851 : iHufnagel et al\ . I2004J : ICohzza et al\ . 12006 J a WM)- They couple en- 



dogenous evolutions of epidemics inside each subpopulation with diffusion 
along the network of connections, which opens interesting new perspectives 
(IColizza et all 12007^ . In this study, we have focused on the issue of arrival 
times of epidemics, as a function of the seed and of the network's character- 
istics. We have proposed a quantity easily computable on any network which 
depends only on the links weights and nodes populations, and which accounts 
for the average arrival time of the spread in each city. This quantity allows to 
sort the various subpopulations according to the order of arrival of each single 
realization of the spread with a very good accuracy. We note that this quan- 
tity is given by a certain weighted distance computed on a directed weighted 
graph, and in particular that this distance selects a shortest weighted path 
between the seed and the various nodes. This result could therefore shed some 
light on the existence and prop erties of "epid e mic pa thways" whose relevant 
role was already suggested by IColizza et al. fl2006J fcl) and which certainly 
deserve further work. In particular, the role and relative importance of the 
multiple paths should be investigated. These predictive tools could play an 
important role in the set-up of containment measures and policies. It would 
also be interesting to adapt th e proposed quantity to more refined or so phisti- 
cated compartmental models (lElveback et all Il976l : IWatts et ali |2005| ) . or to 
generahze it to different s cales such as the urba n scale, where nodes are places 
like home, work or malls (lEubank et all |200J). 
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Fig. 1. Two cities model, distribution of arrival times t^^^ in numerical simulations 
with = 10~^. Inset : distribution of t^^^ with = 10~^ . The continuous line 
is the analytical approximation given by Eq. ([6]). Here and in the following, unless 
specified, N = 10^ and A = 0.1. 



Seed 


{{t),x} 


{{t),Xmp} 


{{t),d} 


{{t),dN/w} 




FRA 


0.856 


0.917 


0.299 


0.311 


0.3 


HKG 


0.874 


0.934 


0.212 


0.275 


0.262 


OTP 


0.866 


0.938 


0.301 


0.304 


0.284 


SJK 


0.883 


0.906 


0.292 


0.287 


0.275 



Table 1 

Kendall's r for the list of cities ranked by arrival order and the list obtained with 
different indicators. Each line corresponds to a different seed: Frankfurt (FRA, de- 
gree 240), Hong-Kong (HKG, degree 109). Bucharest (OTP, degree 56) and Sao 
Jose dos Campos (SJK, Brazil, degree 6). 
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\n(XWw) 

Fig. 2. Slope (A) of (t„) vs n on a line of length 500. The average is taken 
over 500 realizations of the spreading, and the slope is measured at large n (only 
the cities with n > 100 are considered). The error bars represent the variance of A. 
A G [10-2,1] and w G [10, 10^], = 10^. Full line: {A)def Dotted line: {A)Gumbei- 
Dashed line: (A) front- Inset: Correlations Cij = (AjAj) — (Aj)(Aj) versus \ i— j \ 
on a line of length 500, the average is done over 500 realizations. 
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Fig. 3. Distribution of A„ = tn — tn-i for various values of n. For small n, the distri- 
bution is close to a Gumbel and for large n, it evolves to a Gaussian. w/{NX) = 10~^, 
N = 10^ A = 0.1. 
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Fig. 4. Distribution of tn for various values of n. The distribution remains close to 
Gumbel shapes even at large n. The plot on the right is in log-scale on the y-axis 
and the curves are shifted vertically for clarity. w/{NX) = 10~^, = 10^, A = 0.1. 
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Fig. 5. (A-C) Black circles indicate the arrival time distribution on a line at the 
city #7 obtained for a fixed random set of populations {Ni} and weights {wi} 
(3 different sets are used in the 3 graphs). Red crosses correspond to the same 
distribution obtained A) with uniform travel Wi = w, and populations {Ni}; B) 
with uniform populations Ni = N, and weights {wi}; C) with uniform populations 
Ni = N and uniform weights Wi = w. 
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Fig. 6. X{t) vs X for 5 cities connected on a line. Each point correspond to the 
average arrival time in one of the city of the line, for one realization of the weights 
w, averaged over 10^ realizations of the epidemic spread. Results are presented for 
100 different random sets {wi}. 



26 




10 100 

degree of the hub (k) 



Fig. 7. Ratio {ti)/{t\ ) of the average arrival time in a fixed peripheral city i of 
the network shown in the inset and the average arrival time in the Id case, vs the 
degree of the hub, k. Averages are done over 1000 realizations of the spreading, with 
parameter = 10~^. Circles: the seed is the hub; squares: the seed is one of the 
peripheral cities. 
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Fig. 8. Ratio of the average arrival time in B (the seed is A) and the average arrival 
time in the Id case (when only the path ACB is connected), versus w' /w (in the 
simulations wi = W2 = w'l = w'2 = w' ^ and N = N' = Na)- Averages are done 
over 1,000 realizations of the spreading, with w = 10^ and = 10"^. All cities 
have the same population = 10^. Full line: theoretical prediction p7p. 
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Fig. 9. Ratio of the average arrival time in B (the seed is A with degree k) to 
the average arrival time in the Id case (when only the path AlB is connected), 
vs w'/w. Averages are done over 1,000 realizations of the spread, with w = 10'^ 
and = 10~^. All cities have the same population, N = 10^. Numerical results 
(symbols) are compared to the theoretical model given by Eq. (jl9p (full lines) for 
different values of k. 
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Fig. 10. Correlations between the average time {t) of arrival in B and x (B) (black), 
and between {t) and Xmp (B) (red), for the network represented in the inset. The 
different points correspond to different values of k (from 2 to 5), different path 
lengths / (from 2 to 6) and different realizations of the disorder of the links (100 
different sets of {w} for each structure). Averages are done over 1,000 realizations 
of the spreading. 
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Fig. 11. X{t) versus x on the WAN for diseases starting in different cities. Each 
red dot corresponds to a city and averages are done over 1,000 reahzations of the 
spread. Crosses are an average over cities with the same x- The starting seeds are 
Frankfurt (degree 240), Bucharest (degree 56) and Hong-Kong (degree 109). When 
the initial seed is a hub, the average arrival time is larger, especially in the first 
reached cities, due to the effect of the degree. We also show X{t) versus x ^or a SIR 
model seeded in Hong-Kong. 
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Fig. 12. Fraction of couples of nodes correctly ranked as a function of their 
Ax (circles), in each realization of the spread, and cumulative distribution 
(squares) of the values of Ax (i.e., fraction of couples of cities (i,j) with 
^xihj) =1 X U) ~ X (^) l> •^x)- Top: the seed is Hong-Kong, a node of degree 109, 
and the curves of fc for both SI (circles) and SIR (crosses) processes are shown. 
Bottom: the seed has degree 6 (Sao Jose dos Campos, Brazil). 
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Fig. 13. Same quantities as in Fig. [T2l but each graph concerns only couples of 
nodes at the same topological distance from the seed, d. From top left to bottom 
right : d = 2, d = 2>, d = A and d = 5. The seed is Frankfurt (degree 240). 
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Fig. 14. Histogram of Kendall's r for the list of cities sorted according either to the 
arrival time in a given realization or to Different seeds yield different histograms. 
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